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Abstract: Let jj, be a p-dimensional vector, and let Ei and £2 be p x p positive 
definite covariance matrices. On being given random samples of sizes iVi and N2 
from independent multivariate normal populations N p (fj,, £1) and -/V p (/i, £2), re- 
spectively, the Behrens-Fisher problem is to solve the likelihood equations for esti- 
mating the unknown parameters /1, £1, and £2. We shall prove that for Ni,N2 > p 
there are, almost surely, exactly 2p + 1 complex solutions of the likelihood equa- 
tions. For the case in which p = 2, we utilize Monte Carlo simulation to estimate 
the relative frequency with which a typical Behrens-Fisher problem has multiple 
real solutions; we find that multiple real solutions occur infrequently. 

Key words and phrases: Behrens-Fisher problem, Bezout's theorem, maximum like- 
lihood estimation, maximum likelihood degree. 

1 Introduction 

Let /U € W be a p-dimensional vector, and let Si and £2 be p x p positive 
definite (symmetric) matrices. Consider independent multivariate normal popu- 
lations, Np(fj,, Ei) and N p {n, £2), from which we have been given random samples 
Xx, . . . , Xjv x and Y\, . . . , Yn 2 , respectively. On the basis of the given data, the 
famous Behrens-Fisher problem (Behrens (1929), Fisher (1939)) is to estimate 
the parameters n, Si, and £2 by means of the method of maximum likelihood. 

It is well-known that the corresponding system of likelihood equations can- 
not be solved explicitly, and that has led many to propose alternative solutions 
to the Behrens-Fisher problem (Anderson, 2003, p. 187 ff.). More importantly, 
the Behrens-Fisher problem is an early example of a hypothesis testing problem 
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involving exponential families of densities and for which the resulting sufficient 
statistics, when the parameters are restricted to the parameter space determined 
by Ho, fail to be complete (Linnik, 1967). In such a situation, nuisance parame- 
ters exist, and the construction of an exact size-a test is a difficult problem. 

Consequently, the literature on the Behrens-Fisher problem is substantial, 
reflecting the intense interest which the problem has generated since its inception. 
Indeed, the problem has generated an extensive philosophical discussion as well 
as many efforts to derive solutions which are optimal for statistical inference 
(Wallace (1980), Kim and Cohen (1998), Stuart and Ord (1994)). In this paper, 
we determine the number of solutions of the likelihood equations. 

For the case in which p = 1, there are three unknown scalar parameters, 
viz., fi, the common mean, and a\ and a\, the population variances. In this 
case, Sugiura and Gupta (1987) reduced the system of equations to a cubic 
equation in // and deduced that, almost surely, there are three complex solutions; 
they observed also that the likelihood equation tended to have multiple real 
solutions if o~\ and o\ are small in comparison with /i, and otherwise that the 
likelihood equation usually has a unique real solution. Drton (2007) also studied 
the univariate Behrens-Fisher problem and showed, in particular, that if the null 
hypothesis is true then the probability of multiple real solutions tends to zero 
as the sample sizes tend to infinity. We will prove the analogous result for the 
multivariate Behrens-Fisher problem in Theorem 14.11 

In this paper, as in the article of Buot and Richards (2006), we apply results 
from the theory of algebraic geometry to study the solution set of the system of 
likelihood equations for the multivariate Behrens-Fisher problem. Generalizing 
the univariate result described earlier, we shall prove the following result: 

Theorem 1.1. Suppose that N±,N2 > p. Then, almost surely, there are exactly 
2p+l complex solutions of the system of likelihood equations for the multivariate 
Behrens-Fisher problem. In particular, almost surely, there always exists at least 
one real solution. 
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2 Derivation of the likelihood equations 

Denote by X and Y the means of the samples from N p (fi,T,i) and N p (fj,, £2), 
respectively. By standard calculations (cf., Mardia, et al. (1979), p. 142), we 
find that the likelihood equations for estimating fj,, Si and £ 2 are: 

Ni 

t 2 = N^Y J {Y j -Ji){Y j -V)' 
3=1 



(2.1) 



(Nt^ 1 + N 2 ^2 1 )? = N^X + iY 2 £- 1 r. (2.2) 



and 



Some authors have proposed the following iterative algorithm for solving (|2.ip 
and (I2T2D: 



(1) Begin the iteration with initial estimates £j = Si, i = 1, 2, where 

iVj 

(2) Apply (|2.2p to calculate /to, the corresponding estimate of it, in the form 

(3) Use the value of juo obtained in Step (2) to calculate E^i, an updated value 
of E^o, using the formulas 

E Xil = Si + (X- fo)(X - /to)', E 2>1 = S 2 + (Y - po)(F - jEt )', 

which are a consequence of (|2.5|) and (|2.6|) below. 

(4) Return to Step (2) and update jEtj until the sequences Eij and E2J, J = 
1, 2, 3, . . ., converge. 
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We are grateful to Mathias Drton for pointing out that his work in Drton 
and Eichler (2006) implies that this algorithm converges to a saddle point or a 
local (but not necessarily a global) maximum of the likelihood function. If the 
likelihood function were found to be multimodal, a phenomenon which has been 
encountered recently by Drton and Richardson (2004) in a study of seemingly 
unrelated regression models, then any numerical algorithm for solving the sys- 
tem of likelihood equations necessarily must include some information about the 
choice of initial values. 

At first glance, the likelihood equations may appear to be a system of p{p+2) 
equations in p(p+2) variables comprising the p components of fi and the p(p+l)/2 
entries of both Si and £2. However, a closer inspection of (|2.1j) and (|2.2j) reveals 
that if fx is known then Si and £2 are determined completely. We shall show later 
how to eliminate Si and S2 from (|2.2p to obtain a system of p cubic equations 
in the variables jEti, . . . ,fip. 



Proposition 2.1. The likelihood equations \2. 1\) and \2. 2\) for the Behrens- 
Fisher problem are equivalent to 

jViSrHf -ft) + N 2 S^(Y-J1) = Q 

Proof. We apply to the sums in (|2.ip the standard procedure of writing each 
term Xi — fi as Xi — X + X — fi, and similarly for each term Yi — fi. This leads 
to the formulas 

£1 = Si + (X (2.5) 

and 

% = S 2 + {Y -$)(? -V)' (2-6) 

where S\ and S 2 are defined in (|2.3p . By a special case of Woodbury's theorem 
(cf., Muirhead (1982), p. 580, Theorem A5.1) we have, for any nonsingular px p 
matrix M and any column vector v £ MP, 

(M + vv') 1 = M — . 

v ' l + v'M- l v 

Multiplying the latter equation on each side from the right by v and simplifying 
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the result, we obtain 

M- l vv'M~ l v 



(M + vv')~ l v = M^v 



l + v'M- l v 
(1 + v'M- l v)M- 1 v - {M~ l v)(v' M~ l v) 
1 + v'M~ 1 v 

M- l v 



(2.7) 



1 + v'M- 1 v' 

Setting M = S\ and v = X — fl, we obtain 

x-\x-fi) = (s 1 + (x-ji)(x-fiy)-\x-fi) 
i + (x-jiys^(x-jiy 

and, similarly, 

s 2 1 (y-/i) = s 2\Y-ft_ 

l + (Y-i2ys^\Y-Jl) 

On rewriting f)2.2j) as 

Aq£^(X - /2) + iVsE^Y - M) = 0, 
it follows from (j2TTj) and (p^|) that (j2T2|) is equivalent to (j2T4"j) . □ 



3 The maximum likelihood degree of the Behrens- 
Fisher problem 

Following Catanese, et al. (2006) and Ho§ten, et al. (2005) we will call the 
number of complex solutions to the likelihood equations the maximum likelihood 
degree. In this section we shall prove Theorem 11.11 namely, that the maximum 
likelihood (or ML) degree of the Behrens-Fisher problem is 2p + 1. Before pro- 
viding the details of the proof, it is instructive to understand why the theorem 
holds for p = 1 and p = 2. Let us denote by Dxijl) and Dy{p) the denominators 
1 + (X - ty'S^iX - p) and 1 + (Y - p)'S^ l {Y - p), respectively, which appear 
in the likelihood equations (12. 4p . 

Lemma 3.1. Let 

NxDyffiS^iX -p) + N 2 D X (J1)S^(Y - /z) = (3.1) 
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be the system of polynomial equations obtained by clearing denominators in the 
likelihood equations {2.J$ , and suppose that /I is a solution to i3.1\) . Then 
Dx(fi) = if and only if Dy{p) = 0. 

Proof. Suppose that Dx{p) = 0. On multiplying (|3.ip from the left by (X — fi)' 
we obtain N\Dy{fi)(Dx{fi) — 1) = and so we deduce that Dy = 0. Similarly, 
starting with the assumption that Dy = 0, we deduce that Dx = 0. □ 

We remark that, because Dx(p) and Dy(p) are strictly positive for any real 
ju, the system of equations (|2.4p and (13. ip are equivalent when determining real 
solutions only. However, in the calculation of complex solutions, the likelihood 
equations (|2.4p are not equivalent to (|3.ip since it is possible that the denomina- 
tors are zero for complex fi. 

Let J be the ideal defined by the equations (|3.ip and let I = (Dj(/i), Dy(fl)) 
be the ideal of zeros common to the first and second denominators. Then we need 
to compute and count the solutions to J : I. For the case in which p = 1 there is a 
single univariate cubic polynomial in (]3.ip which, generically, has three complex 
roots. Since two generic univariate polynomials (in this case, Dx{fi) and Dy{fi)) 
have no common roots then the ideal I has, in general, no solutions. Hence we 
conclude that J : I has exactly three solutions for the case in which p = 1. 

We now consider the case in which p = 2. Since two quadrics in two variables 
have, generically, four complex roots then there are four generic solutions to I. 
Similarly, since two cubics in two variables have generically nine complex roots 
then there are nine generic solutions to J. Therefore J : I has five complex roots 
for the case in which p = 2. 

Unfortunately, this counting argument fails even for p = 3. In this case, we 
have two quadrics in three variables, so there are infinitely many solutions to / 
and hence also to J. Yet, J : I still has finitely many solutions. Theorem 11.11 
follows from the following result of Catanese, et al. (2006). 

Theorem 3.2. (Catanese, et al, 2006) Let fi, . . . , f n be polynomials of degrees 
b\, . . . , b n , respectively, in the variables x\, . . . , x^; let u\, . . . ,u n be integers; let 
f = Z" 1 • • • f% n ; and consider the critical equations 

10/ _ 13/ = ... = 19/ = Q 

/ dxx f dx 2 f dx n 
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of log / = X^r=i u i 1°§ fi- V th e number of complex solutions to this system of 
equations is finite then that number is less than or equal to the coefficient of z d 
in the generating function 



Equality holds if the coefficients of the polynomials fi are sufficiently generic. 

Before we proceed, there are a few points that need clarification in Theorem 
13.21 First of all, given the integers b\,...,b n there exists a fixed polynomial 
G = G% _ _ _ 5 n in the coefficients of n polynomials in d variables with degrees 
61, . . . , b n so that we call fi, ■ ■ ■ , f n generic if G(/i, . . . , f n ) 7^ 0. Furthermore, 
when /1 , . . . , f n is generic in this sense the number of complex solutions to the 
critical equations is given by the formula in the statement of the theorem. In 
other words, genericity already implies the finiteness of the number of complex 
solutions. This follows from Theorem 5 in Catanese, et al. (2006) leading up to 
the proof of Theorem 13.21 

Before providing the proof of Theorem 11.11 we show that the coefficients of 
Dxiji) an d Dy(p) are generic for almost all data X\, . . . , X^ 1 and Yi, . . . , Y/v 2 . 
First we need the following result which has a standard proof in the literature 
(for instance based on the argument on page 76 in Anderson (2003)). We present 
our own proof. 

Lemma 3.3. Suppose that N + 1 > p. Then given any p x p positive definite 
matrix S there exist X x , . . . , X N+1 € W such that S = Y,i=i l ( X i ~ x )( x i ~ X Y ■ 

Proof. We can assume, without loss of generality that X = 0. Now let Xi = 
(Xn, . . . , Xip)' for i = 1, . . . , N and let X^+i = — J2iLi x i- Given a posi- 
tive definite matrix S, there exists a nonsingular symmetric matrix U such that 
USU' = A where A is a diagonal matrix with diagonal entries Aj > 0, i = 1, . . . ,p. 
Hence it is enough to prove the result for diagonal matrices A. The required iden- 
tity A = X^S 1 x i X i gives rise to p(p + l)/2 polynomial equations, namely, 




(l-b 1 z){l-b 2 z)---{l-b n z)' 




N N 



i=l j=l 
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k = 1, . . . ,p, and 

N N N 

Xu(X\j + Xfcj) + X 2 i(X 2 j + Xjy) + h Xm(XNj + -Xfcj) = 0, 

fc=i fe=i fc=i 

1 < i < i < P- We claim that there exists at least one real solution to the 
above system where Xy = for j = 1, . . . ,p and £ = j + 1, . . . , N. It is easy 
to check that X u = ^/4± with X lk = X 2k = ■■■ = X k _ lik = J k( ^+i) and 
X kk = -kJ ^^jy for k = 2, . . . , p give such a solution. □ 

Theorem 3.4. For generic data X\, . . . , Xx 1 and Yi, . . . , Y^ 2 the denominators 
Dx and Dy are generic. 

Proof. Dx and Dy are quadratic forms in p variables. In the light of our remarks 
after Theorem 13.21 there exists a fixed polynomial G in the coefficients of two 
quadratic forms in p variables such that Dx and Dy are generic if G(Dx, Dy) ^ 
0. We need to show that this condition holds for generic data. As in the proof of 
Lemma [3.3[ the entries of S\ and S 2 are polynomials in the data. The same lemma 
implies that the polynomial maps defined from the data spaces RP xAfl and W xNi 
are surjective onto the cone of semidefinite matrices in Rp(p +1 )/ 2 . Therefore there 
exist data vectors X\, . . . , Xx x and Yi, . . . , Yx 2 such that G(Dx, Dy) ^ 0. If the 
statement in the theorem is not true, then there exists a Zariski open subset 
U C W xNl x RP xN * such that for all (X 1 ,...,X m : Yi,...,Y Na ) € U we 
have G(Dx , Dy) = 0. But this means that G is identically zero, and this is a 
contradiction. □ 

Proof of Theorem, li.il Denoting by .LQu, £1, E2) the likelihood function for the 
Behrens-Fisher problem, then it is well-known that 

L(£,Ei,E 2 ) = (2vre)-( iVl + iV2 ^/ 2 |S 1 r JVl / 2 |S 2 r 7V2 / 2 . 
By (|2.5p and (|2.6p . we have [Ej| = \Si\ -Dxijl) an d |£ 2 | = 1 5*2 1 ■ Dy (p) . Therefore 
L(/2,Ei,S 2 ) 

= (2vre)-( iVl+JV2 ^ 2 \Si\- N ^ 2 \S 2 \- N ^ 2 {Dx(P)Y Nl ' 2 {Dy{p)y N2/2 . 
It now is clear that, to find the maximum value of L, we need to minimize 

+ nys^(x - n)) N ^ 2 (1 + (y - pys 2 -\Y - n)) N2/2 . (3.2) 
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Equivalently, we may minimize the logarithm of this expression, and since the 
critical equations of the logarithm of (13, 2p are precisely the likelihood equations 
in Proposition ^, 11 then Theorem 13 . 2 1 implies that the maximum likelihood degree 
of the Behrens-Fisher problem is equal to the coefficient of z p in the power series 
expansion of the rational function 



provided that the data X\, . . . , and Yi, . . . , Yjv 2 > and hence Dx and Dy are 
generic. By expanding this rational function in a power series in z, we find that 
this coefficient equals 



and an elementary calculation shows that this sum equals 2p + 1. 

Step 4 in Algorithm 7 in Ho§ten, et al. (2005) and the theory of Grobner 
bases imply that all 2p + 1 complex solutions can be obtained from the roots 
of a univariate polynomial of degree 2p + 1. Since Dx(fi) and Dyifi) have real 
coefficients, this univariate polynomial also has real coefficients. In particular, 
since roots occur in complex conjugate pairs then at least one root is real. □ 

Remark 3.5. We note that our arguments which led to the derivation of the ML 
degree of the Behrens-Fisher problem also apply to the more general problem of 
multivariate analysis of variance (MANOVA). Suppose that we have independent 
multivariate normal populations N p (fi, Si), . . . , N p (/j,, £fc+i) and that, on the ba- 
sis of random samples from each population, we wish to derive the maximum like- 
lihood estimators of the parameters ji and Si, ... , £fc+i. By arguments similar 
to those in Section [H we obtain analogous likelihood equations as in Proposition 
12.11 where now there are k + 1 rational summands in each of the p equations. It 
then follows from Theorem 13.21 that the ML degree for the MANOVA problem is 



(l-2z) 2 ' 



i+j=p 





(3.3) 



By writing this result in the form 



p 



p\ ^j + k 



) 



d{k, P ) = i+Y J {-±) p - j y 



3=1 
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we find that d(k,p) is odd; therefore, there always exists a real solution to the 
system of likelihood equations. 

We note that d(k,p) can be evaluated using methods from the calculation of 
combinatorial sums, as follows: First, we write 



Inserting this formula in the sum in (|3.3p and interchanging derivatives and 
summation, we obtain 



In particular, d(l,p) = 2p + l, the ML degree of the Behrens-Fisher problem, and 
d(2,p) = 2p(p + 1) + 1. The general formula for d(k,p) is interesting even in the 
case p = 1, for it yields the ML degree of the one-dimensional {k + l)-population 
MANOVA problem to be 2k + 1. Further, by substituting t = (1 + u)/2 in (pT4")) . 
we recognize the outcome as Rodrigues' formula (Szego, 1939, p. 66) for a Jacobi 
polynomial pjf fc , and we obtain d(k,p) = pjf fc,0 ^(3), p> k. 

4 Simulations and a large sample size result 

Having determined the number of solutions of the system of likelihood equations 
(|3.ip it is natural to seek the number of real solutions, for it is those solutions 
which are of interest in statistical inference. Not surprisingly, it appears to be 
difficult to determine an algebraic expression for the number of real solutions 
of the system; indeed, this is also the case for the general theory of systems of 
polynomial equations. 

To study the real solutions of the system (13.11) . we considered the case in 
which p = 2, presenting empirical evidence that multiple solutions occur rarely if 
the model is correctly specified. In each simulation run, we first used a random 
number generator to generate sample sizes A^i and N2, and a mean vector [i. 
We next generated lower triangular matrices T\ and T2 with positive diagonal 





(3.4) 
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entries, after which we set = T^T^, k = 1,2. Finally, we simulated a random 
sample of vectors Z±,... , Zn x from A^(0, 12), and then we set Xj = T\Zj + 
j = 1, . . . , N\. It follows from standard distribution theory that X\, . . . , Xn 1 
constitutes a simulated sample from the bivariate normal population A^(M)Si). 
In a similar manner, we simulated an independent random sample Y±, . . . , 
from N 2 (n, £2)- 

The solutions of the resulting likelihood equations (|3.ip were computed nu- 
merically using PHCpack (Verschelde, 1999), a software package which implements 
polyhedral homotopy continuation methods for solving systems of polynomial 
equations. The results of our simulations show that multiple solutions can occur. 
For example, for N% = 11, N2 = 5, and the summary statistics 

% _ ( -1-5516 \ S ( ' 3 " 8 ~°- 1026 ^ 

y -9.4713 ) ' 1 y -0.1026 0.2378 J ' 

^ / -1.9175 \ ~ / 0.4193 0.0792 \ 

y -10.4805 / ' ' 1 \ 0.0792 0.0334 J ' 

the real solutions for \i are 

I -1.3570 \ ( -1.2478 \ ^ t -1.4451 \ 
\ -10.2957 J \ -9.9902 J \ -9.6333 j 

This example seems, however, to be a rare exception. Indeed, we found that the 
bivariate Behrens-Fisher likelihood equations (|3.1|) had one real solution in about 
99.5% of simulations, three real solutions in about 0.5% of simulations, and we 
found no instances in which the equations had five real solutions. However, it 
is possible that (|3.ip has five real solutions when the data is generated from a 
"wild" distribution and not from the corresponding multivariate distributions. 
For instance, for N\ = 15, N2 = 28, and 



X = _ , Si 




Y = , 5 : 



49.3619 -45.0547 

-45.0547 42.4495 

52.8534 19.8380 \ 

19.8380 9.0472 / 
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the real solutions for fj, are 

/ 3.9822 \ / -3.7286 \ / -2.4192 \ / 2.0437 \ / 1.0089 \ 
y 1.0443 J ' y 3.2906 J ' \ 4.6925 J ' \ 5.8993 J ' \ 8.2001 J ' 

To test for distinctions between the case of small and large samples in the 
bivariate case, we performed simulations in which N\ and N2 were randomly 
generated (uniform distribution) between 3 and 15. The outcomes are given as 
follows, with percentages rounded-off to two decimal places: 



Please place Table 1 here 



As noted above, none of these simulation resulted in five real solutions. 

In the case of larger samples, our simulations resulted in the following out- 
comes: 



Please place Table 2 here 



Here again, no simulation resulted in five real solutions. (In both cases, the 
population mean fi is randomly generated from a uniform distribution on the 
subspace [—20,20] x [—20,20], and the population covariance matrices Si and 
£2 are randomly generated in the manner described above, with positive diagonal 
entries whose values are no greater than 10.) 

In summary, there seems to be little chance that a randomly generated, two- 
dimensional Behr ens-Fisher problem will have three or more real solutions, and 
there is a high chance that it will have a unique real solution. The following 
supports the second conclusion for large sample sizes. 

Theorem 4.1. Suppose that the random samples X±, . . . , Xn 1 and Y±,..., Yn 2 

are drawn from independent normal populations N p (fi, Si) and N p (p,, S2), re- 
spectively. As Ni,N% — > 00 the likelihood equations \2.J$ for the Behrens- Fisher 
problem has a unique real root with probability one. 



Proof. If X = Y then it follows from (13. 2j) that the unique real solution of the 
likelihood equations is /2 = X = Y . Without loss of generality we can assume 
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that X = Y = 0, and with this the likelihood equations are 

N ^p + N *Sp = 0. (4.1) 
1 + ju'S^ju 1 + p'S Y l p 

We argue that ju = is a solution of multiplicity one for the system obtained by 
clearing denominators in (|4.ip . Let / be the ideal in C[p.±, . . . ,p p ] generated by 
these p equations. The multiplicity of p, = is the length of the artinian module 

C[/Xj, ■ ■ ■ ^p]^!,.. .,/J.p) 

I-C[fix,.. ■ 

over the local ring C[/ii, . . . , /%>]( M1) ... )/Lt }■ I is generated by p polynomials given 
by (A^ 1 S x 1 + A^25 y 1 )^ + A^i(^ , 5 y 1 /I)^ 1 ^ + iV2(/I'5 x 1 /i)5 y 1 /I = 0. Each of these 
polynomials consists of a linear term and a cubic term. With probability one the 
rank of NiS^ 1 + /^.Sy 1 over C is p, and hence we can assume that / is generated 
by p polynomials of the form fa + gi where gi has degree three. This implies that 
the initial ideal of I in the local ring C[/*i, . . . , Pp]^,...,^ ) with respect to the local 
term order anti-graded revlex as on page 152 of Cox, et al. (1998) is (p±, . . . , fj, p ). 
By Corollary 4.5 of Cox, et al. (1998), we conclude that the length of the above 
module and hence the multiplicity of /I = is one. Now as Ni,N^ — > oo, by the 
Law of Large Numbers, X and Y converge to p, and S\ and S2 converge to Sx 
and Sy- Since /2 = is the unique real solution to (|4.ip with multiplicity one, 
and by the continuity of solutions to the general likelihood equations (|2.4p . we 
conclude that with probability one (|2,4p has a unique real solution. □ 
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Table 1: Simulations with 3 < N 1 , N 2 < 15 



Number of solutions 


Frequency 


Percentage 


1 

3 


4450 
32 


99.29% 
0.71% 



Table 2: Simulations with 15 < N x , N 2 < 60 



Number of solutions 


Frequency 


Percentage 


1 


4404 


99.46% 


3 


24 


0.54% 
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